{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single-Qubit Noisy Simulator\n",
    "\n",
    "\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    " This tutorial describes how to use Quanlse's single-qubit noisy simulator to define the decoherence noise and amplitude noise and simulate the gate operation. The outline of this tutorial is as follows:\n",
    " \n",
    "- Introduction\n",
    "- Preparation\n",
    "- Define simulator parameters and pulses\n",
    "- $\\pi$ pulse calibration using Rabi oscillation \n",
    "- Measure relaxation time $T_1$ \n",
    "- Measure decoherence time $T_2$ by Ramsey experiment\n",
    "- Single-qubit noisy simulator at the gate level\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "At the pulse level, a single-qubit's control process requires applying pulses to the qubit. In this tutorial, the single-qubit noisy simulator enables the simulation of the dynamical evolution of the control based on the control pulse and noise parameters input by the user.\n",
    "\n",
    "The system Hamiltonian of a three-level superconducting qubit in the rotating frame after RWA (Rotating Wave Approximation) can be written as \\[1\\]:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm sys}(t) = \\hat{H}_{\\rm anharm} + \\hat{H}_{\\rm noise}(t) + \\hat{H}_{\\rm ctrl}(t),\n",
    "$$\n",
    "\n",
    "where the anharmonicity term $\\hat{H}_{\\rm anharm}$ is related to qubit's anharmonicity parameter $\\alpha$:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm anharm} = \\frac{\\alpha}{2} \\hat{a}^\\dagger \\hat{a}^\\dagger \\hat{a}\\hat{a},\n",
    "$$\n",
    "\n",
    "where $\\hat{a}^\\dagger$ and $\\hat{a}$ are the three-level qubit creation operator and annihilation operator, respectively. In the single-qubit noisy simulator, we introduce two kinds of noise: decoherence noise due to the environment and amplitude noise due to fluctuations of the control pulse waveform $\\hat{H}_{\\rm amp}$ \\[1\\]:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm amp}(t) = \\epsilon \\hat{H}_{\\rm ctrl},\n",
    "$$\n",
    "\n",
    "where the probability distribution of the amplitude noise parameter $\\epsilon$ follows Gaussian distribution:\n",
    "\n",
    "$$\n",
    "P(\\epsilon) = \\frac{1}{\\sqrt{ 2 \\pi \\sigma^2 }}e^{-\\frac{ \\epsilon^2} {2 \\sigma^2} }.\n",
    "$$\n",
    "\n",
    "The parameters $T_1$ and $T_2$ related to decoherence describe the decay rate of non-diagonal elements of the density matrix and the population decay rate on the excited state. We can use Bloch-Redfield density matrix $\\rho_{BR}$ to illustrate the decoherence process of an initialized qubit $|\\psi\\rangle = \\alpha |0\\rangle + \\beta |1\\rangle$ under open-system evolution ($\\delta \\omega$ is the difference between the frequency of the control pulse and the frequency of the qubit) \\[2\\]:\n",
    "\n",
    "$$\n",
    "\\rho_{BR} =\n",
    "\\begin{pmatrix}\n",
    "1 + (|\\alpha|^2-1)e^{-t/T_1} & \\alpha\\beta^* e^{i\\delta \\omega t}e^{-t/T_2} \\\\\n",
    "\\alpha^* \\beta e^{-i\\delta\\omega t} e^{-t/T_2} & |\\beta|^2 e^{-t/T_1}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "As it shows, the decay rate of the population on the excited state is related to $T_1$, and the decay rate of the non-diagonal elements is related to $T_2$.\n",
    "\n",
    "Therefore, we use the coefficients $\\sigma$, $T_1$, and $T_2$ to characterize these two kinds of noise, which means the user can change the three parameters to study the noises of the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "After you have successfully installed Quanlse, you could run the Quanlse program below following this tutorial. To run this particular tutorial, you would need to import the following packages from Quanlse and other commonly-used Python libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import 1-qubit noisy simulator at the pulse level\n",
    "from Quanlse.Superconduct.Simulator import PulseModel\n",
    "from Quanlse.Superconduct.Simulator.PulseSim1Q import pulseSim1Q\n",
    "from Quanlse.QOperator import driveX, driveZ\n",
    "from Quanlse.QWaveform import gaussian, square\n",
    "from Quanlse.Utils.Functions import basis, dagger, expect, project\n",
    "from Quanlse.QOperation.FixedGate import H, X, Y, Z\n",
    "from Quanlse.QOperation.RotationGate import RX, RY, RZ\n",
    "from Quanlse.Utils.Bloch import plotBloch, rho2Coordinate\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "\n",
    "# Import tool for analysis\n",
    "from scipy.optimize import curve_fit\n",
    "from matplotlib import pyplot as plt\n",
    "from math import pi\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the Quanlse Cloud Service, we need to acquire a token from http://quantum-hub.baidu.com to get access to the cloud. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Define class and set the token for cloud service\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define simulator parameters and pulses\n",
    "\n",
    "We first define the required parameters of the simulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubitNum = 1  # The number of qubits\n",
    "dt = 0.2  # The sampling period, in nanoseconds\n",
    "level = 3  # The energy level\n",
    "anharm = -0.3472 * (2 * pi)  # The anharmonicity of the qubit, in 2 * pi * GHz\n",
    "wq = 4.9  * (2 * pi)  # The qubit frequency, in 2 * pi * GHz\n",
    "\n",
    "# Define the noise of the simulator\n",
    "ampSigma = 0.02  # amplitude (over-rotation) error\n",
    "t1 = 3000  # qubit relaxation time, in nanoseconds\n",
    "t2 = 800  # qubit dephasing time, in nanoseconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an object of class `PulseModel`. The object is a physics model defined by the parameters above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubitAnharm = {0: anharm}\n",
    "qubitFreq  = {0: wq}\n",
    "qubitT1 = {0: t1}\n",
    "qubitT2 = {0: t2}\n",
    "\n",
    "model = PulseModel(subSysNum=qubitNum, sysLevel=level, dt=dt, ampSigma=ampSigma,\n",
    "                   T1=qubitT1, T2=qubitT2, qubitFreq=qubitFreq, qubitAnharm=qubitAnharm)\n",
    "\n",
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have constructed the physics model of the simulator. Then we can use Rabi oscillation to find the amplitudes of $\\pi$ and $\\pi / 2$ pulse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\pi$ pulse calibration using Rabi oscillation\n",
    "\n",
    "In this section, we use Rabi oscillation to calibrate $\\pi$ pulse. Namely, we fix the duration time of the pulse, then vary the amplitude of the driving pulse and record the population of state $|1\\rangle$. This process can be done by the single-qubit simulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the amplitudes of the pulse\n",
    "ampList = np.linspace(0, 0.4, 100)\n",
    "\n",
    "# Create a jobList object\n",
    "rabiJob = ham.createJobList()\n",
    "\n",
    "# Define the shape of the gaussian waveform\n",
    "tg = 60\n",
    "tau = tg / 2\n",
    "sigma = tg / 8\n",
    "\n",
    "# Append each job of different pulse amplitudes to jobList\n",
    "for amp in ampList:\n",
    "    wave = gaussian(t=tg, a=amp, tau=tau, sigma=sigma)\n",
    "    job = ham.createJob()\n",
    "    job.appendWave(operators=driveX, onSubSys=0, waves=wave)\n",
    "    job = model.getSimJob(job)\n",
    "    rabiJob.addJob(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the `jobList` of Rabi oscillation, we call the module `runHamiltonian()` with input initial state $|\\psi\\rangle = |0\\rangle$ and the joblist. We set the parameter `isOpen=True` to simulate the evolution in the open system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the initial state for the simulation\n",
    "stateInit = basis(level, 0) \n",
    "    \n",
    "# Run the simulation of the open system evolution\n",
    "result = runHamiltonian(ham=ham, state0=stateInit, jobList=rabiJob, isOpen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we plot the population on state $|1\\rangle$ as a function of the driving pulse amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the projector\n",
    "prj = basis(level, 1) @ dagger(basis(level, 1))\n",
    "\n",
    "popList = []\n",
    "\n",
    "# Compute the population for each job\n",
    "for res in result:\n",
    "    rho = res['state']\n",
    "    popList.append(expect(prj, rho))  # Compute the population of |1>\n",
    "\n",
    "plt.plot(ampList, popList, '.')\n",
    "plt.xlabel('Amplitudes (a.u.)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.title('Rabi oscillation')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the corresponding cosine function by fitting these points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the function to be fitted\n",
    "def fit(x, a, b, c, d):\n",
    "    return a * np.cos(b * x + c) + d\n",
    "\n",
    "# Fit the curve\n",
    "paraFit, _ = curve_fit(fit, ampList, popList, [-0.5, 2 * np.pi / 0.3, 0, 0.5])\n",
    "def yFit(x):\n",
    "    return fit(x, paraFit[0], paraFit[1], paraFit[2], paraFit[3])\n",
    "y = [yFit(x) for x in ampList]\n",
    "\n",
    "# Plot the fitted curve\n",
    "plt.plot(ampList, y)\n",
    "plt.xlabel('Amplitudes (a.u.)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.ylim(-0.05, 1)\n",
    "plt.title('Rabi oscillation')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it shows, the population on $|1\\rangle$ oscillates periodically with the increase of pulse amplitude. Using the fitted function, we can find the corresponding amplitudes when the population of $|1\\rangle$ is 0.5 and 1, which are exactly the amplitudes of $\\pi / 2$ and $\\pi$ pulse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ampList = np.linspace(0, 0.4, 5000)\n",
    "\n",
    "piAmp = []\n",
    "halfPiAmp = []\n",
    "for amp in ampList:\n",
    "    if abs(yFit(amp) - 0.5) < 1e-3:\n",
    "        halfPiAmp.append(amp)\n",
    "    if abs(yFit(amp) - 0.98) < 1e-3:\n",
    "        piAmp.append(amp)\n",
    "\n",
    "# find the corresponding amplitudes\n",
    "x90 = min(halfPiAmp)\n",
    "x180 = min(piAmp)\n",
    "\n",
    "# Print the results\n",
    "print(f'The amplitudes of pi/2 pulse: {x90}')\n",
    "print(f'The amplitudes of pi pulse: {x180}')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measure relaxation time $T_1$ \n",
    "\n",
    "To measure $T_1$ relaxation time, we first implement a pi pulse to maximize the population on the excited state. Then we record the evolution of the population over time to obtain the relaxation time $T_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the idling time\n",
    "tIdle = 2 * t1\n",
    "\n",
    "# Define the wave to flip the qubit state\n",
    "wave = gaussian(t=tg, a=x180, tau=tau, sigma=sigma)\n",
    "\n",
    "# Initialize a job\n",
    "job = ham.createJob()\n",
    "\n",
    "# Firstly, apply a X gate to flip the qubit\n",
    "job.appendWave(operators=driveX, onSubSys=0, waves=wave)\n",
    "\n",
    "# Then simulate the evolution during the idling time\n",
    "job.appendWave(operators=driveX, onSubSys=0, waves=square(tIdle, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the task and the initial state, we run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the initial state for the simulation\n",
    "stateInit = basis(level, 0) \n",
    "    \n",
    "# Run the simulation of the open system evolution\n",
    "result = runHamiltonian(ham=ham, state0=stateInit, job=job, isOpen=True)\n",
    "\n",
    "# Define the projector |1><1|\n",
    "prj = basis(level, 1) @ dagger(basis(level, 1))\n",
    "\n",
    "popList = []\n",
    "\n",
    "# Calculate the population of |1> during the evolution\n",
    "for rho in result[0]['evolution_history']:\n",
    "    popList.append(expect(prj, rho))\n",
    "\n",
    "# Get the maximum time of the job\n",
    "maxTime, _ = job.computeMaxTime()\n",
    "\n",
    "tList = np.linspace(0, maxTime, len(popList))\n",
    "\n",
    "# Plot the time-evolution poplulation for simulation and prediction \n",
    "plt.plot(tList, popList, '-', label='simulation')\n",
    "tIdleList = np.linspace(tg, tIdle, 20)\n",
    "plt.plot(tIdleList, np.exp(-1. / t1 * np.array(tIdleList - tg)), label='prediction')\n",
    "plt.xlabel('Time (ns)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.title(r'$T_1$ measurement')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure above is the dynamical evolution of population of the excited state $|1\\rangle$. It can be seen that within tens of nanoseconds from the beginning of the experiment, the applied pulse of the $X$ gate drives the qubit from the ground state to the excited state. Then in the following idling time, the decay rate of the excited state's population is consistent with the theoretical value. At the time $t=T_1$, the population of $|1\\rangle$ decays to $1/e$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measure decoherence time $T_2$  by Ramsey experiment\n",
    "\n",
    "To measure decoherence time $T_2$, we first implement a $\\pi / 2$ pulse. After a delay time $t_{\\rm idle}$, we implement another $\\pi / 2$ pulse again and obtain the population of $|1\\rangle$. In this process, we can observe oscillation over the delay time due to the detuning between the frequency of drive pulse and qubit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the maximum idling time\n",
    "maxTime = t2\n",
    "\n",
    "# Define the detuning \n",
    "detuning = 2 * pi * 8. / maxTime\n",
    "\n",
    "# Define the job for Ramsey experiment\n",
    "tList = np.linspace(0, maxTime, 200)\n",
    "ramseyJob = ham.createJobList()\n",
    "\n",
    "for t in tList:\n",
    "    job = ham.createJob()\n",
    "    job.appendWave(driveX, 0, gaussian(t=tg, a=x90, tau=tau, sigma=sigma), compact=False)  # pi/2 pulse\n",
    "    job.appendWave(driveZ, 0, waves=square(t, detuning), compact=False)  # simulate the rotation due to the detuning\n",
    "    job.appendWave(driveX, 0, gaussian(t=tg, a=-x90, tau=tau, sigma=sigma), compact=False)  # pi/2 pulse\n",
    "    ramseyJob.addJob(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the list of jobs and initial state, we run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the simulation starting with initial state |0>.\n",
    "stateInit = basis(level, 0)\n",
    "result = runHamiltonian(ham, stateInit, jobList=ramseyJob, isOpen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the population of $|1\\rangle$ over the delayed time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popList = []\n",
    "\n",
    "prj = basis(level, 1) @ dagger(basis(level, 1))\n",
    "\n",
    "for res in result:\n",
    "    rho = res['state']\n",
    "    popList.append(expect(prj, rho))\n",
    "\n",
    "plt.plot(tList, popList, '.b', label='simulation')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we fit the curve with a defined fitting function, and estimate $T_2$ according to the parameters of the fitting function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the fitting function\n",
    "def fitRamsey(x, a, b):\n",
    "    return - np.cos(a * x) * np.exp(- b * x) * 0.5 + 0.5\n",
    "\n",
    "paraFit, _ = curve_fit(fitRamsey, tList, popList, [detuning, 0.])\n",
    "\n",
    "def yFit(x):\n",
    "    return fitRamsey(x, paraFit[0], paraFit[1])\n",
    "\n",
    "# Plot the result\n",
    "plt.plot(tList, popList, '.b', label='simulation')\n",
    "plt.plot(tList, yFit(tList), label='fit')\n",
    "plt.plot(tList, np.exp(- (1 / t2 + 1/(2 * t1)) * tList) * 0.5 + 0.5, label='prediction')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Print the estimated T_2 time\n",
    "print(f'The T2-dephasing time is approximately {1 / (paraFit[1] - 1 / (2 * t1))} ns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single-qubit noisy simulator at the gate level\n",
    "\n",
    "Besides, we can also use the simulator directly at the gate level. First, we call a prepared single-qubit object of `PulseSim1Q()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pulseSim1Q(dt=0.2)\n",
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define a set of quantum gates. Here, we firstly define a set of single-qubit gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H(model.Q[0])\n",
    "X(model.Q[0])\n",
    "Y(model.Q[0])\n",
    "Z(model.Q[0])\n",
    "H(model.Q[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use method `model.schedule()` to generate the corresponding pulse sequence, and call `runHamiltonian` to simulate the evolution of single-qubit with the presence of decoherence. We can visualize the process via module `plotBloch()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = model.schedule()\n",
    "\n",
    "res = runHamiltonian(ham, state0=basis(3, 0), job=job, isOpen=True)\n",
    "\n",
    "history = res[0]['evolution_history']\n",
    "\n",
    "posList = []\n",
    "\n",
    "for rho in history:\n",
    "    rho2d = project(rho, 1, 3, 2) / np.trace(rho)\n",
    "    posList.append(rho2Coordinate(rho2d))\n",
    "\n",
    "plotBloch(posList, save=True, mode='animate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it shows, with the decoherence, the Bloch vector of qubit would not stay on the surface of the sphere, indicating the pure initial state has evolved to a mixed state, and the non-diagonal elements of the density matrix decayed. The system lost information due to the interaction with the environment. Therefore, the relatively small $T_1$ and $T_2$ will reduce the performance of deep quantum circuits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This tutorial describes how to simulate a superconducting single-qubit process considering partial noise using Quanlse and visualize the results. Users can click on this link [tutorial-single-qubit-noisy-simulator.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/EN/tutorial-single-qubit-noisy-simulator.ipynb) to jump to the corresponding GitHub page for this Jupyter Notebook documentation to get the relevant code, try the different parameter values for further exploring the function of the Quanlse Simulator module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\n",
    "\\[1\\] [Carvalho, Andre RR, et al. \"Error-robust quantum logic optimization using a cloud quantum computer interface.\" *arXiv preprint arXiv:2010.08057* (2020).](https://arxiv.org/abs/2010.08057)\n",
    "\n",
    "\\[2\\] [Krantz, Philip, et al. \"A quantum engineer's guide to superconducting qubits.\" *Applied Physics Reviews* 6.2 (2019): 021318.](https://aip.scitation.org/doi/abs/10.1063/1.5089550)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
